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ABSTRACT 

In this paper we present a supervised hyperspectral image 
segmentation algorithm based on a convex formulation of 
a marginal maximum a posteriori segmentation with hid¬ 
den fields and structure tensor regularization: Segmentation 
via the Constraint Split Augmented Lagrangian Shrinkage 
by Structure Tensor Regularization (SegSALSA-STR). This 
formulation avoids the generally discrete nature of segmen¬ 
tation problems and the inherent NP-hardness of the integer 
optimization associated. 

We extend the Segmentation via the Constraint Split Aug¬ 
mented Lagrangian Shrinkage (SegSALSA) algorithm (H by 
generalizing the vectorial total variation prior using a struc¬ 
ture tensor prior constructed from a patch-based Jacobian (2). 
The resulting algorithm is convex, time-efficient and highly 
parallelizable. This shows the potential of combining hid¬ 
den fields with convex optimization through the inclusion of 
different regularizers. The SegSALSA-STR algorithm is val¬ 
idated in the segmentation of real hyperspectral images. 

Index Terms — Image segmentation, hidden fields, struc¬ 
ture tensor regularization, Constrained Split Augmented La¬ 
grangian Shrinkage Algorithm (SALSA) 

1. INTRODUCTION 

Supervised image segmentation is fundamental in a large 
number of hyperspectral image applications m. Image seg¬ 
mentation aims to partition the image domain such that pixels 
belonging to the same partition element share similar proper¬ 
ties, either by raw pixel values or more complex features. Due 
to its ill-posed nature ( e.g . depending on the goal, multiple 
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similarity criteria for grouping pixels in the same partition 
element can be used), image segmentation is often performed 
with help of regularization that promote or penalize different 
behaviors, e.g. the use of a Markov Random Field (4l MRF) 
to promote smoothness of the labeling. 

As partitions are naturally represented by images of in¬ 
tegers, maximum a posteriori (MAP) segmentations become 
integer optimization problems. SegSALSA cm sidesteps 
the discrete nature of the segmentation problems by using 
a hidden field to drive the segmentation and marginalizing 
on the hidden field m combined with a vectorial total vari- 
ation (VTV) prior 0 0 . The segmentation is inferred by 
computing the marginal maximum a posteriori (MMAP) es¬ 
timate of the hidden field, which is a convex problem solved 
by the Split Augmented Lagrangian Shrinkage (SALSA) al¬ 
gorithm EJ. 

In this paper, we extend the formulation of SegSALSA 
to include a generalization of VTV prior based on structure 
tensor priors 0, thus deriving a more general algorithm for 
hyperspectral image segmentation. The structure tensor prior 
is constructed from a patch-based Jacobian. The regulariza¬ 
tion via the structure tensor prior is achieved via a Schatten 
norm regularization, akin to ED. 

The paper is organized as follows. Section [2] describes 
the SegSALSA algorithm, introducing the hidden fields and 
the marginal MAP formulation. Section [3] presents the struc¬ 
ture tensor regularization and describes the construction of 
the structure tensor prior from the patch-based Jacobian. Sec¬ 
tion [4] formulates our optimization problem and presents the 
SegSALSA-STR algorithm. Section [5] illustrates the validity 
of the algorithm on synthetic and real hyperspectral data, and 
Section [6] concludes the paper. 

2. BACKGROUND 

Let x G R dxn denote a n-pixel hyperspectral image with d 
spectral bands, where G R d denotes the feature vector cor¬ 
responding to the image pixel i, S = {1, ..., n} the set that 



indexes the image pixels, C = { 1 ,..., K} the set of possible 
K labels, and y E C n a labelling of the image. 


Maximum a posteriori segmentation The SegS ALSA for¬ 
mulation adopts a Bayesian perspective; the MAP segmenta¬ 
tion associated with the labeling y is obtained as 

y = arg max p(y|x) = arg max p(x|y)p(y), ( 1 ) 

ye L N 


where p(y |x) and p(x|y) denote the posterior probability and 
the observation model respectively, and p( y) the prior on the 
labeling. Under assumption of conditional independence of 
the observation model, we can expand the observation model 
p(x|y) = UiesPMVi)- The optimization associated MAP 
formulation Q is an integer optimization problem: apart 
from a small number of exceptions, the use of contextual 
priors p( y), such as MRF priors makes 0 a combinatorial 
problem. 

Hidden fields and marginal MAP To mitigate the diffi¬ 
culties associated with the integer optimization problem of 
the MAP, we move from the discrete formulation by intro¬ 
ducing a hidden field ( 6 | and marginalizing on the discrete 
labels. We represent the hidden field z as a K x n ma¬ 
trix containing a collection of hidden vectors z x E R K for 
i E S. The joint probability of the labeling y and the field 
z is p( y,z) = p(y|z)p(z), with assumption of conditional 
independence of p(y |z) = YliesP(yi\ z i)- The joint prob¬ 
ability of the features x, the labeling y and the field z is 
p(x, y, z) = p(x|y)p(y|z)p(z). We can now marginalize on 
the discrete labels, 

f>( x , z ) = n { ^2 P( x i>2/iMy*> z i)M z )> 

i£S 

and obtain a marginal MAP (MMAP) estimate of the hidden 
field, 


ZMMAP = arg max p(x, z) 

= arg min — lnp(xlz) — lnp(z), (2) 

z£M Kxn 

which is no longer a discrete optimization problem. 

Link between class labels and hidden fields As a model 
for the conditional probabilities p(yijzi) for the labels given 
the field, we adopt 


p(yi = k\zi) = [z i] k , 

for i e S and k E C. This link between the probabilities of 
the class labels and the hidden fields imposes two constraints 
on the hidden field and, consequently, on the optimization 0 : 
the field must be nonnegative, and each hidden vector z x for 


i E S must sum to one. We can now formulate the optimiza¬ 
tion ([ 2 | as 

zmmap = arg min - lnp(x|z) - lnp(z), (3) 

zeM Kxn 

subject to: z > 0 , l^z = l n . 


3. STRUCTURE TENSOR REGULARIZATION 

We extend the SegSALSA algorithm by replacing the VTV 
prior by a generalization of VTV based on structure tensor 
regularization 0 . 


Patch-based Jacobian Following closely the notation 
in m, we define the patch-based Jacobian of the hidden 
field as 


[Mf = 


(PiD h z)f 
(PiD v z)f 


(■ P L D h z)f- 

(P l D v z)T\ ’ 


( 4 ) 


where [Jz\i denotes the components of the patch-based Jaco¬ 
bian on the ith pixel of the hyperspectral image, correspond¬ 
ing to a ( KL ) x 2 matrix. D^ and D v are the horizontal 
and vertical difference operators, and Pj is a weighted shift 
operator corresponding to the patch, with each operator be¬ 
ing applied equally to the entire field D h ,D v , Pj : R Kxn 
R Kxn . Assuming the patch to be a (2M+1) x (2M+1) rect¬ 
angular patch, with L = (2 M + l ) 2 pixels, Pj corresponds 
to the jth possible shift within the patch (from the L possible 
shifts) weighted by a Gaussian centered on the center of the 
patch and with a bandwidth 7 . 


Discrete structure tensor From the patch-based Jacobian 
of the hidden field 0, the structure tensor Sl is defined as 

[S L z]i = [Jz]f [Jz]i, (5) 

a 2 x 2 matrix for the Ah pixel of the hyperspectral image. 

The minimization of the eigenvalues of the structure ten¬ 
sor Sl in 0 leads to the penalization of variations of the 
field among the pixels in patch. As there is an intrinsic con¬ 
nection between the eigenvalues of the structure tensor 0 
(A+, A_) and the singular values of the patch-based Jacobian 
(a/^+ 5 we can minimize the singular values instead. 

Let || [ Jz]i ||s p denote the Schatten p norm of the patch-based 
Jacobian 

||[J Z y Sp = \\a{[Jz] i )\\ p , 

where a([Jz]i) represent the singular values of [Jz]*. The 
discrete structure tensor prior can be constructed through the 
minimization of the singular values of the patch-based Jaco¬ 
bian 

-lnp(z) = Xj2\\[Jz\i\\s p +c te . (6) 

It should be noted that for 1 x 1 patches ( L = 1), the min¬ 
imization of the Schatten norm of the structure tensor is 
equivalent to the minimization of the VTV, leading to the 
SegSALSA formulation in [DEI- 




4. OPTIMIZATION ALGORITHM 

With the SegSALSA general formulation described in Section 
[ 2 ] and armed with the new generalized total variation prior 
from the structure tensor described in Section [3] we can now 
describe our algorithm which we name SegSALSA-STR. 


Salsa methodology The optimization (\T\ is solved follow¬ 
ing the SALSA methodology (9), an instance of the alternat¬ 
ing direction method of multipliers designed to optimize sum 
of an arbitrary number of convex terms. Let d denote the 
scaled Lagrange multipliers; solving ClD reduces to the fol¬ 
lowing iterative decoupled problem 


Problem formulation Combining the MMAP formulation 
in ([3]) with the prior described in ([ 6 ]), we can formulate our 
MMAP problem as follows: 

zmmap =argmin Y] - ln(pf z*) + A Y] || [Jz];|| s (7) 

zeM KxnZ — z ' 

subject to: z >0 l^z = 1 ^, 

where p; = \pfc\yi = 1).. .p(x|^ = K )] T , and pfz* > 0 
in the feasible set. As the Hessian of — ln(pf z^) is semidef- 
inite positive, and ||[Jz]i||s p is a composition of norms, the 
optimization 0 is convex. The solution z M map is computed 
by the SALSA methodology 0. However, we cannot apply 
the SALSA methodology directly to 0 . To do so we rewrite 
the problem as a sum of convex functions with linear con¬ 
straints. 


Rewriting the optimization problem We rewrite the opti¬ 
mization problem ([7]) as 


4 


min 


rfXn 


3 = 1 


( 8 ) 


where Hj, for j = 1,..., 4, denote linear operators, and g^, 
for j = 1,..., 4 denote closed, proper, and convex functions. 
We define the linear operators Hj as: 


Hr ml, H 2 = J , H 3 = I, H 4 = I, (9) 

where I : 'R Kxn —>> 'R Kxn denotes the identity operator and 
J : R Kxn M 2 LKxn denotes an operator stacking the 
patch-based Jacobians defined in ([?]). We define the closed, 
proper, and convex functions g 3 as: 

9i(0 = X]- ln (pfO. 92(0 = ||£||s p , (10) 

93(0 = i +(0, 94(0 = iiiixO- 

The functions ?’ + and i\ are indicator functions for the sets 
xn and l n , respectively. The introduction of a variable 
splitting 

Ujf = Hjz, 

for j = 1,..., 4, and with u E M( 3 + 2L ) Kxn 9 allows us to re¬ 
formulate the optimization ^ into a constrained formulation 

4 

min ^ gj (uj ), subject to u = Gz, (11) 

U ’ Z 3=1 

with G the corresponding to the stacking of the linear opera¬ 
tors Hj, for j = 1,..., 4. 


z fc+1 = argmin ||Gz - u k - d k \\ 2 F , ( 12 ) 

z 

4 

u fe+1 = argmin^^(u J ) + |||Gz fe -u-d fc ||2 r , (13) 

3 = 1 

d k+i = d fc _ ( Gz fc+1 _ u fe+i). ( 14 ) 


Solving the iterative decoupled problem The quadratic 
problem ( [ 12 ] ) is solved by computing independent cyclic 
convolutions on each image of z, with time complexity 
O(LKnlogn). The problem (T3| ) can be decoupled for 
each of the convex functions gj, amounting to computing 
the Moreau proximity operator (MPO) CD for each of the 
functions. 

For j = 1,3,4, involved in the respective MPO is respec¬ 
tive the MPOs, all with O(Kri) complexity are as follows: g\ 
is the root of a polynomial using a closed form solution; g 3 is 
the projection onto the positive orthant; and g 4 is the projec¬ 
tion onto the probability simplex. For j = 2, the MPO is the 
solution of the following problem 


u. 


fc+i 


= arg min 

u 2 £R( 2LK ) x n 


E 


IlIuAlk 


F \\H 2 z k -\i 2 -d%\\ 2 F . 


For p = 1, which we adopt in this paper, this can be solved 
by the soft thresholding of the singular values of H 2 z k — d k . 


Complexity, parallelization and stopping criterion The 

time complexity of solving the problem ( f]~3j ) is dominated 
by the computation of n single value decompositions of ma¬ 
trices (LK) x 2, which amounts to a time complexity of 
0{n(LK ) 3 ) El. The complexity of SegSALSA-STR is 

0(LKn\ogn + n(LK ) 3 ), 

being highly parallelizable: LK decoupled fast Fourier trans¬ 
forms followed by n decoupled SVDs. As a stopping crite¬ 
rion, we impose that both the primal and dual residuals are 
smaller than a given threshold. It was observed that a fixed 
number of iterations of the order of 100 provides excellent 
results. 


5. RESULTS 

We validate our algorithm in the segmentation of a subscene 
of the AVIRIS Indian Pine scene (Fig. [lj. The subscene con¬ 
sists of a 145 x 145 pixel subsection with 200 spectral bands 
and contains 16 different classes. 15 Samples per class are 



Fig. 1 . AVIRIS Indian Pines scene. SegSALSA-STR segmentation results for multiple patch sizes, starting from the same MLR segmen¬ 
tation. Top row — (a) ground truth, ( b ) MLR segmentation - 56.86% accuracy, (c) segmentation for 5 x 5 patches - 78.60% accuracy, (d) 
hidden field for grass-trees class (light blue) for segmentation (c), (e) hidden field for soybean-mintill class (dark yellow) for segmentation 
(c). Bottom row — (/) segmentation for 1 x 1 patches (corresponding to SegSALSA with VTV priors) - 74.33% accuracy, ( g ) segmentation 
for 3 x 3 patches — 77.53% accuracy, (/i) segmentation for 7x7 patches — 78.64% accuracy, (i) segmentation for 9 x 9 patches — 78.33%, 
(j) segmentation for 11 x 11 patches - 78.01% accuracy. 


used as training set. The classes are modeled by a multinomial 
logistic regression (MLR) and the MLR weights are learnt us¬ 
ing the LORSAL 113J. The parameter A corresponding to the 
relative weight 0 of the structure tensor prior is set to 2 , a 
compromise that allows good classification performance with 
a lower number of iterations of the algorithm. 

We illustrate the joint effect of the size of the patch and 
of the bandwidth 7 of the Gaussian weighting for the con¬ 
struction of the patch-based Jacobian 0 , by comparing the 
segmentation results of a real hyperspectral image (Fig. [TJ 
for varying sizes of the patches. The band with 7 is set as the 
distance in pixels from the border of the patch to the center 
(e.g. 7 = 2 for 5 x 5 patches), with 7 = 1 for 1 x 1 patches. 
As structure tensor regularization is a generalization of the 
VTV prior, in the case of 1 x 1 patches, SegSALSA-STR cor¬ 
responds to SegSALSA with the VTV prior U0. 

The analysis of the segmentation results by SegSALSA- 
STR shows the capabilities of the structure tensor prior in ob¬ 
taining smooth segmentation boundaries. The hidden fields 
corresponding to two different classes (Figs, d, e) illustrate 
the preservation of detail of the segmentation (Fig. c). The 
effect of increasing the patch size is evident (bottom row) as 
the segmentation performance for higher patch sizes is higher 
than the segmentation performance for lower patch sizes. 

The use of structure tensor priors with a patch size larger 
than lxl (Figs, c, g, h , i ), allows an increase of segmenta¬ 
tion performance when compared to 1 x 1 patches (Fig. /). 
This shows that extending the SegSALSA algorithm by using 


a generalization of the VTV prior based on structure tensor 
regularization allows for improved segmentation performance 
at the cost of a higher computational burden. 

6. CONCLUDING REMARKS 

In this paper we extended the SegSALSA algorithm to use a 
generalized total variation prior based on structure tensor reg¬ 
ularization — SegSALSA-STR. This is a more general for¬ 
mulation of the SegSALSA algorithm which extends the con¬ 
cept of VTV from the pixel level to the patch level. The shift 
from a discrete to a continuous formulation paves the way for 
a class of relevant approaches in remote sensing through the 
inclusion of different priors. The algorithm was validated in 
simulated and real hyperspectral data. 
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